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Quantum annealing is emerging as a promising near-term quantum computing approach to solv¬ 
ing combinatorial optimization problems. A solver for the job-shop scheduling problem that makes 
use of a quantum annealer is presented in detail. Inspired by methods used for constraint satisfac¬ 
tion problem (CSP) formulations, we first define the makespan-minimization problem as a series of 
decision instances before casting each instance into a time-indexed quadratic unconstrained binary 
optimization. Several pre-processing and graph-embedding strategies are employed to compile opti¬ 
mally parametrized families of problems for scheduling instances on the D-Wave Systems’ Vesuvius 
quantum annealer (D-Wave Two). Problem simplifications and partitioning algorithms, including 
variable pruning, are discussed and the results from the processor are compared against classical 
global-optimum solvers. 


I. I. INTRODUCTION 

The commercialization and independent benchmark¬ 
ing m of quantum annealers based on superconduct¬ 
ing qubits has sparked a surge of interest for near-term 
practical applications of quantum analog computation in 
the optimization research community. Many of the early 
proposals for running useful problems arising in space 
science [S] have been adapted and have seen small-scale 
testing on the D-Wave Two processor [B]. The best proce¬ 
dure for comparison of quantum analog performance with 
traditional digital methods is still under debate 13 HI IS] 
and remains mostly speculative due to the limited num¬ 
ber of qubits on the currently available hardware. While 
waiting for the technology to scale up to more significant 
sizes, there is an increasing interest in the identification 
of small problems which are nevertheless computation¬ 
ally challenging and useful. One approach in this direc¬ 
tion has been pursued in [!J|, and consisted in identifying 
parametrized ensembles of random instances of opera¬ 
tional planning problems of increasing sizes that can be 
shown to be on the verge of a solvable-unsolvable phase 
transition. This condition should be sufficient to observe 
an asymptotic exponential scaling of runtimes, even for 
instances of relatively small size, potentially testable on 
current- or next-generation D-Wave hardware. An em¬ 
pirical takeaway from |B] (validated also by experimen¬ 
tal results in [10j [TT]) was that the established program¬ 
ming and program running techniques for quantum an¬ 
nealers seem to be particularly amenable to scheduling 
problems, allowing for an efficient mapping and good 
performance compared to other applied problem classes 
like automated navigation and Bayesian-network struc¬ 
ture learning [Si- 

Motivated by these first results, and with the inten¬ 
tion to challenge current technologies on hard problems 
of practical value, we herein formulate a quantum an¬ 
nealing version of the job-shop scheduling problem (JSP). 
The JSP is essentially a general paradigmatic constraint 
satisfaction problem (CSP) framework for the problem 


of optimizing the allocation of resources required for the 
execution of sequences of operations with constraints on 
location and time. We provide compilation and running 
strategies for this problem using original and traditional 
techniques for parametrizing ensembles of instances. Re¬ 
sults from the D-Wave Two are compared with classical 
exact solvers. The JSP has earned a reputation for be¬ 
ing especially intractable, a claim supported by the fact 
that the best general-purpose solvers (CPLEX, Gurobi 
Optimizer, SCIP) struggle with instances as small as 10 
machines and 10 jobs (10 x 10) [JTT| . Indeed, some known 
20 x 15 instances often used for benchmarking still have 
not been solved to optimality even by the best special- 
purpose solvers El. and 20 x 20 instances are typically 
completely intractable. We note that this early work con¬ 
stitutes a wide-ranging survey of possible techniques and 
research directions and leave a more in-depth exploration 
of these topics for future work. 

A. Problem definition and conventions 

Typically the JSP consists of a set of jobs J = 
{j i,..., jjv} that must be scheduled on a set of machines 
M = {mi,..., mu/}- Each job consists of a sequence of 
operations that must be performed in a predefined order 

jn = {O n \ —> O n 2 —>■■■—> OnL n }. 

Job j„ is assumed to have L n operations. Each opera¬ 
tion O n j has an integer execution time p n j (a value of 
zero is allowed) and has to be executed by an assigned 
machine m g . £ AJ, where q n j is the index of the as¬ 
signed machine. There can only be one operation run¬ 
ning on any given machine at any given point in time and 
each operation of a job needs to complete before the fol¬ 
lowing one can start. The usual objective is to schedule 
all operations in a valid sequence while minimizing the 
makespan (i.e., the completion time of the last running 
job), although other objective functions can be used. In 
what follows, we will denote with T the minimum possi¬ 
ble makespan associated with a given JSP instance. 
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As defined above, the JSP variant we consider is de¬ 
noted 3M\p nj £ [Pmin,... ,_Pmax]|C' max in the well-known 
a|/9|7 notation, where p min and p max are the smallest and 
largest execution times allowed, respectively. In this no¬ 
tation, J M stands for job-shop type on M machines, and 
Cmax means we are optimizing the makespan. 

For notational convenience, we enumerate the opera¬ 
tions in a lexicographical order in such a way that 


jl 

= {Oi —►•••—»• 

°k J, 

h 

= {Ofc 1+ 1 —>■••• 

“A Ok 2 }, 

j N 

= {Ok N _ 1 +1 -A 

• • • Ok N } 


Given the running index over all operations i £ 
{1,..., /cat}, we let qi be the index of the machine m g! re¬ 
sponsible for executing operation Oi. We define I. m to be 
the set of indices of all of the operations that have to be 
executed on machine m m , i.e., I m = {i : qi = m}. The 
execution time of operation Oi is now simply denoted pi. 

A priori, a job can use the same machine more than 
once, or use only a fraction of the M available machines. 
For benchmarking purposes, it is customary to restrict a 
study to the problems of a specific family. In this work, 
we define a ratio 9 that specifies the fraction of the total 
number of machines that is used by each job, assuming 
no repetition when 9 < 1. For example, a ratio of 0.5 
means that each job uses only 0.5 M distinct machines. 


B. Quantum annealing formulation 

In this work, we seek a suitable formulation of the JSP 
for a quantum annealing optimizer (such as the D-Wave 
Two). The optimizer is best described as an oracle that 
solves an Ising problem with a given probability m- 
This Ising problem is equivalent to a quadratic uncon¬ 
strained binary optimization (QUBO) problem [TD]. The 
binary polynomial associated with a QUBO problem can 
be depicted as a graph, with nodes representing variables 
and values attached to nodes and edges representing lin¬ 
ear and quadratic terms, respectively. The QUBO solver 
can similarly be represented as a graph where nodes rep¬ 
resents qubits and edges represent the allowed connec¬ 
tivity. The optimizer is expected to find the global min¬ 
imum with some probability which itself depends on the 
problem and the device’s parameters. The device is not 
an ideal oracle: its limitations, with regard to precision, 
connectivity, and number of variables, must be consid¬ 
ered to achieve the best possible results. As is customary, 
we rely on the classical procedure known as embedding 
to adapt the connectivity of the solver to the problem 
at hand. This procedure is described in a number of 
quantum annealing papers lain. During this proce¬ 
dure, two or more variables can be forced to take on 
the same value by including additional constraints in the 
model. In the underlying Ising model, this is achieved 
by introducing a large ferromagnetic (negative) coupling 


Jf between two spins. The embedding process modifies 
the QUBO problem accordingly and one should not con¬ 
fuse the logical QUBO problem value, which depends on 
the QUBO problem and the state considered, with the 
Ising problem energy seen by the optimizer (which addi¬ 
tionally depends on the extra constraints and the solver’s 
parameters, such as Jp). 

We distinguish between the optimization version of the 
JSP, in which we seek a valid schedule with a mini¬ 
mal makespan, and the decision version, which is lim¬ 
ited to validating whether or not a solution exists with a 
makespan smaller than or equal to a user-specified time- 
span T. We focus exclusively on the decision version 
and later describe how to implement a full optimization 
version based on a binary search. We note that the deci¬ 
sion formulation where jobs are constrained to fixed time 
windows is sometimes referred in the literature as the 
job-shop CSP formulation |16j Hz], and our study will 
refer to those instances where the jobs share a common 
deadline T. 


II. II. QUBO PROBLEM FORMULATION 

While there are several ways the JSP can be formu¬ 
lated, such as the rank-based formulation [18] or the dis¬ 
junctive formulation m , our formulation is based on a 
straightforward time-indexed representation particularly 
amenable to quantum annealers (a comparative study of 
mappings for planning and scheduling problems can be 
found in my We assign a set of binary variables for each 
operation, corresponding to the various possible discrete 
starting times the operation can have: 

{ 1 : operation Oi starts at time t, 

0 : otherwise. ' ' 

Here t is bounded from above by the timespan T, which 
represents the maximum time we allow for the jobs to 
complete. The timespan itself is bounded from above by 
the total work of the problem, that is, the sum of the 
execution times of all operations. 


A. Constraints 


We account for the various constraints by adding 
penalty terms to the QUBO problem. For example, an 
operation must start once and only once, leading to the 
constraint and associated penalty function 



= 1 for each i 


ED 


'i,t 



(3) 


There can only be one job running on each machine at 
any given point in time, which expressed as quadratic 
constraints yields 


E 




= 0 for each m, 


(4) 


3 



operation 1 

operation 2 

operation 3 

jl 

mi, p = 2 

m 2 , p = 1 

m 3 , p = 1 

J 2 

m 3 , p = 2 

mi, p = 1 

m 2 , p = 2 

J3 

m 2 , p = 1 

m 1( p = 1 

m 3 , p = 2 



Machine 1 Machine 2 Machine 3 

d) 

h 
t 2 
U 

*4 

t 5 
*6 
tj 

Machine 1 Machine 2 Machine 3 

FIG. 1: a) Table representation of an example 3x3 instance 
whose execution times have been randomly selected to be ei¬ 
ther 1 or 2 time units, b) Pictorial view of the QUBO map¬ 
ping of the above example for Ht= e- Green, purple, and 
cyan edges refer respectively to hi, / 12 , and /13 quadratic cou¬ 
pling terms (Eqs. [7}]9|). Each circle represents a bit with its 
i, t index as in Eq. 72] c) The same QUBO problem as in 
(b) after the variable pruning procedure detailed in the sec¬ 
tion on QUBO formulation refinements. Isolated qubits are 
bits with fixed assignments that can be eliminated from the 
final QUBO problem, d) The same QUBO problem as in 
(b) for Ht= 7 . Previously displayed edges in the above figure 
are omitted. Red edges/circles represent the variations with 
respect to Ht= 6 - Yellow stars indicate the bits which are 
penalized with local fields for timespan discrimination. 



wlrere R m = A m U B m and 

Am — { [%, t, k, t ) . (t, /t) €: Im X Im: 

i ^ k, 0 < t, t' < T, 0 < t' — t < pi], 

B m — {[%, t, k, t ) . (i, k) €: Im X Im, 

i < k,t' = t,pi > 0 ,pj > 0}. 


The set A rn is defined so that the constraint forbids oper¬ 
ation Oj from starting at t' if there is another operation 
Oi still running, which happens if Oi started at time t 
and t' — t is less than pi. The set B m is defined so that 
two jobs cannot start at the same time, unless at least 
one of them has an execution time equal to zero. Finally, 
the order of the operations within a job are enforced with 


E 

k n — i<i<.k n 
t+pi>t r 




for each n, 


(5) 


which counts the number of precedence violations be¬ 
tween consecutive operations only. 

The resulting classical objective function (Hamilto¬ 
nian) is given by 


H t {x) = ijhi{x) + ah 2 {x ) + /3h 3 (x), 


where 


friO) = E 


( 


E 




k n — 1 <i<k n 

\ t+pi>t' 


h 2 (a;) = E E x i,tXk,f 

h 3{x) = [e^ - ^ ' 


( 6 ) 


(7) 

( 8 ) 
(9) 


and the penalty constants 77 , a, and /3 are required to 
be larger than 0 to ensure that unfeasible solutions do 
not have a lower energy than the ground state (s). As ex¬ 
pected for a decision problem, we note that the minimum 
of Ht is 0 and it is only reached if a schedule satisfies all 
of the constraints. The index of Ht explicitly shows the 
dependence of the Hamiltonian on the timespan T , which 
affects the number of variables involved. Figure [l}b il¬ 
lustrates the QUBO problem mapping for Ht =6 for a 
particular 3x3 example (Figure [lj a). 


B. Simple variable pruning 

Figure [l]b also reveals that a significant number of the 
NMT binary variables required for the mapping can be 
pruned by applying simple restrictions on the time index 
t (whose computation is polynomial as the system size in¬ 
creases and therefore trivial here). Namely, we can define 
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an effective release time for each operation correspond¬ 
ing to the sum of the execution times of the preceding 
operations in the same job. A similar upper bound cor¬ 
responding to the timespan minus all of the execution 
times of the subsequent operations of the same job can 
be set. The bits corresponding to these invalid starting 
times can be eliminated from the QUBO problem alto¬ 
gether since any valid solution would require them to be 
strictly zero. This simplification eliminates an estimated 
number of variables equal to NM (M ( p ) — 1), where (p) 
represents the average execution time of the operations. 
This result can be generalized to consider the previously 
defined ratio 9, such that the total number of variables 
required after this simple QUBO problem pre-processing 
is 6NM[T - 0M(p) + 1], 

III. III. QUBO FORMULATION REFINEMENTS 

Although the above formulation proves sufficient for 
running JSPs on the D-Wave machine, we explore a few 
potential refinements. The first pushes the limit of simple 
variable pruning by considering more advanced criteria 
for reducing the possible execution window of each task. 
The second refinement proposes a compromise between 
the decision version of the JSP and a full optimization 
version. 


A. Window shaving 

In the time-index formalism, reducing the execution 
windows of operations (i.e., shaving) [20], or in the dis¬ 
junctive approach, adjusting the heads and tails of opera¬ 
tions E 31211, or more generally, by applying constraints 
propagation techniques (e.g. [12), together constitute 

the basis for a number of classical approaches to solving 
the JSP. Shaving is sometimes used as a pre-processing 
step or as a way to obtain a lower bound on the makespan 
before applying other methods. The interest from our 
perspective is to showcase how such classical techniques 
remain relevant, without straying from our quantum an¬ 
nealing approach, when applied to the problem of prun¬ 
ing as many variables as possible. This enables larger 
problems to be considered and improves the success rate 
of embeddability in general (see Figure [3]), without sig¬ 
nificantly affecting the order of magnitude of the overall 
time to solution in the asymptotic regime. Further im¬ 
mediate advantages of reducing the required number of 
qubits become apparent during the compilation of JSP 
instances for the D-Wave device due to the associated 
embedding overhead that depends directly on the num¬ 
ber of variables. The shaving process is typically handled 
by a classical algorithm whose worst-case complexity re¬ 
mains polynomial. While this does not negatively impact 
the fundamental complexity of solving JSP instances, for 
pragmatic benchmarking the execution time needs to be 
taken into account and added to the quantum annealing 


runtime to properly report the time to solution of the 
whole algorithm. 

Different elimination rules can be applied. We focus 
herein on the iterated Carlier and Pinson (ICP) proce¬ 
dure m reviewed in the appendix with worst-case com¬ 
plexity given by 0{N 2 M 2 T\og(N)). Instead of looking 
at the one-job sub-problems and their constraints to elim¬ 
inate variables, as we did for the simple pruning, we look 
at the one-machine sub-problems and their associated 
constraints to further prune variables. An example of 
the resulting QUBO problem is presented in Figure [I[c. 


B. Timespan discrimination 

We explore a method of extracting more information 
regarding the actual optimal makespan of a problem 
within a single call to the solver by breaking the de¬ 
generacy of the ground states and spreading them over 
some finite energy scale, distinguishing the energy of valid 
schedules on the basis of their makespan. Taken to the 
extreme, this approach would amount to solving the full 
optimization problem. We find that the resulting QUBO 
problem is poorly suited to a solver with limited preci¬ 
sion, so a balance must be struck between extra informa¬ 
tion and the precision requirement. A systematic study 
of how best to balance the amount of information ob¬ 
tained versus the extra cost will be the subject of future 
work. 

We propose to add a number of linear terms, or local 
fields, to the QUBO problem to slightly penalize valid 
solutions with larger makespans. We do this by adding 
a cost to the last operation of each job, that is, the set 
{Ofcj,... ,Ok N }- At the same time, we require that the 
new range of energy over which the feasible solutions are 
spread stays within the minimum logical QUBO prob¬ 
lem’s gap given by A E = minjry, a, /3}. If the solver’s 
precision can accomodate K distinguishable energy bins, 
then makespans within [T — K, T] can be immediately 
identified from their energy values. The procedure is il¬ 
lustrated in Figure [l}d and some implications are dis¬ 
cussed in the appendix. 


IV. IV. ENSEMBLE 
PRE-CHARACTERIZATION AND 
COMPILATION 

We now turn to a few important elements of our com¬ 
putational strategy for solving JSP instances. We first 
show how a careful pre-characterization of classes of ran¬ 
dom JSP instances, representative of the problems to be 
run on the quantum optimizer, provides very useful in¬ 
formation regarding the shape of the search space for T. 
We then describe how instances are compiled to run on 
the actual hardware. 
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A. Makespan Estimation 

In Figure [2] we show the distributions of the opti¬ 
mal makespans T for different ensembles of instances 
parametrized by their size N = M, by the possible val¬ 
ues of task durations V p = {p m i n) • • • ,Tmax}j and by the 
ratio 9 < 1 of the number of machines used by each 
job. Instances are generated randomly by selecting OM 
distinct machines for each job and assigning an execution 
time to each operation randomly. For each set of parame¬ 
ters, we can compute solutions with a classical exhaustive 
solver in order to identify the median of the distribution 
(T)(N,V pi 9) as well as the other quantiles. These could 
also be inferred from previously solved instances with 
the proposed annealing solver. The resulting informa¬ 
tion can be used to guide the binary search required to 
solve the optimization problem. Figure [2] indicates that 
a normal distribution is an adequate approximation, so 
we need only to estimate its average (T) and variance 
a 2 . Interestingly, from the characterization of the fam¬ 
ilies of instances up to N = 10 we find that, at least 
in the region explored, the average minimum makespan 
(T) is proportional to the average execution time of a 
job ( p)9N , where (p) is the mean of V v . This linear 
ansatz allows for the extrapolation of approximate re¬ 
source requirements for classes of problems which have 
not yet been pre-characterized, and it constitutes an ed¬ 
ucated guess for classes of problems which cannot be pre¬ 
characterized due to their difficulty or size. The accu¬ 
racy of these functional forms was verified by computing 
the relative error of the prediction versus the fit of the 
makespan distribution of each parametrized family up to 
N = M = 9 and p max = 20 using 200 instances to com¬ 
pute the makespan histogram. The prediction for (T) 
results are consistently at least 95% accurate, while the 
one for a has at worst a 30% error margin, a very ap¬ 
proximate but sufficient model for the current purpose of 
guiding the binary search. 


B. Compilation 

The graph-minor embedding technique (abbreviated 
simply “embedding”) represents the de facto method of 
recasting the Ising problems to a form compatible with 
the layout of the annealer’s architecture [MJ J25], which 
for the D-Wave Two is a Chimera graph pQ. Formally, 
we seek an isomorphism between the problem’s QUBO 
graph and a graph minor of the solver. This procedure 
has become a standard in solving applied problems using 
quantum annealing El HU and can be thought of as the 
analogue of compilation in a digital computer program¬ 
ming framework during which variables are assigned to 
hardware registers and memory locations. This process 
is covered in more details in the appendix. An example of 
embedding for a 5 x 5 JSP instance with 9=1 and T = 7 
is shown in Figure [3]-a, where the 72 logical variables of 
the QUBO problem are embedded using 257 qubits of the 
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FIG. 2: a) Normalized histograms of optimal makespans T 
for parametrized families of JSP instances with N = M, 
V p on the y-axis, 9 = 1 (yellow), and 9 = 0.5 (purple). 
The distributions are histograms of occurrences for 1000 ran¬ 
dom instances, fitted with a Gaussian function of mean (T). 
We note that the width of the distributions increases as the 
range of the execution times V p increases, for fixed (p). The 
mean and the variance are well fitted respectively by (T) = 
A T A Pm in T / 1; N Pmax and <7 — Cq T C a \ T) T A CT p tn i ti + B a Pmax, 
with At = 0.67, Bt = 0.82, ao = 0.7, A a = —0.03, 
B a = 0.43, and C a = 0.003. 


Chimera graph. Finding the optimal tiling that uses the 
fewest qubits is NP-hard [25], and the standard approach 
is to employ heuristic algorithms m- In general, for the 
embedding of time-indexed mixed-integer programming 
QUBO problems of size N into a graph of degree k, one 
should expect a quadratic overhead in the number of bi¬ 
nary variables on the order of aN 2 , with a < (fc — 2) _1 
depending on the embedding algorithm and the hardware 
connectivity m- This quadratic scaling is apparent in 
Figure [3}b where we report on the compilation attempts 
using the algorithm in m- Results are presented for the 
D-Wave chip installed at NASA Ames at the time of this 
study, for a larger chip with the same size of Chimera 
block and connectivity pattern (like the latest chip cur¬ 
rently being manufactured by D-Wave Systems), and for 
a speculative yet-larger chip where the Chimera block 
is twice as large. We deem a JSP instance embeddable 
when the respective Ht=t is embeddable, so the decrease 
in probability of embedding with increasing system size is 
closely related to the shift and spreading of the optimal 
makespan distributions for ensembles of increasing size 
(see Figure 2]). What we observe is that, with the avail¬ 
able algorithms, the current architecture admits embed¬ 
ded JSP instances whose total execution time NM9(p) 
is around 20 time units, while near-future (we estimate 2 
years) D-Wave chip architectures could potentially dou¬ 
ble that. As noted in similar studies (e.g., (6]) , graph 
connectivity has a much more dramatic impact on em- 
beddability than qubit count. 
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FIG. 3: a) Example of an embedded JSP instance on NASA’s 
D-Wave Two chip. Each chain of qubits is colored to represent 
a logical binary variable determined by the embedding. For 
clarity, active connections between the qubits are not shown, 
b) Embedding probability as a function of N = M for 9 = 1 
(similar results are observed for 6 = 0.5). Solid lines refer to 
V p = [1,1] and dashed lines refer to V p = [0, 2]. 1000 random 
instances have been generated for each point, and a cutoff 
of 2 minutes has been set for the heuristic algorithm to find 
a valid topological embedding. Results for different sizes of 
Chimera are shown, c) Optimal parameter-setting analysis 
for the ensembles of JSP instances we studied. Each point 
corresponds to the number of qubits and the optimal Jf (see 
main text) of a random instance, and each color represents a 
parametrized ensemble (green: 3x3, purple: 4x4, yellow: 
5x5, blue: 6x6; darker colors represent ensembles with V v = 
[1,1] as opposed to lighter colors which indicate V p = [0, 2]). 
Distributions on the right of scatter plots represent Gaussian 
fits of the histogram of the optimal Jf for each ensemble. 
Runtime results are averaged over an ungauged run and 4 
additional runs with random gauges [25]. 


Once the topological aspect of embedding has been 
solved, we set the ferromagnetic interactions needed to 
adapt the connectivity of the solver to the problem at 
hand. For the purpose of this work, this should be re¬ 


garded as a technicality necessary to tune the perfor¬ 
mance of the experimental analog device and we include 
the results for completeness. Introductory details about 
the procedure can be found in ia m- In Figure [3}c 
we show a characterization of the ensemble of JSP in¬ 
stances (parametrized by N, M, 9, and V p , as described 
at the beginning of this section). We present the best 
ferromagnetic couplings found by runs on the D-Wave 
machine under the simplification of a uniform ferromag¬ 
netic coupling by solving the embedded problems with 
values of Jp from 0.4 to 1.8 in relative energy units of the 
largest coupling of the original Ising problem. The run 
parameters used to determine the best Jf are the same 
as we report in the following sections, and the problem 
sets tested correspond to Hamiltonians whose timespan 
is equal to the sought makespan Ht=t- This parameter¬ 
setting approach is similar to the one followed in [B] for 
operational planning problems, where the instance en¬ 
sembles were classified by problem size before compila¬ 
tion. What emerges from this preliminary analysis is 
that each parametrized ensemble can be associated to a 
distribution of optimal Jp that can be quite wide, espe¬ 
cially for the ensembles with p m i n = 0 and large p max . 
This spread might discourage the use of the mean value 
of such a distribution as a predictor of the best Jp to 
use for the embedding of new untested instances. How¬ 
ever, the results from this predictor appear to be better 
than the more intuitive prediction obtained by correlat¬ 
ing the number of qubits after compilation with the op¬ 
timal Jp. This means that for the D-Wave machine to 
achieve optimal performance on structured problems, it 
seems to be beneficial to use the information contained in 
the structure of the logical problem to determine the best 
parameters. We note that this “offline” parameter-setting 
could be used in combination with “online” performance 
estimation methods such as the ones described in [28] in 
order to reach the best possible instance-specific Jp with 
a series of quick experimental runs. The application of 
these techniques, together with the testing of alternative 
offline predictors, will be the subject of future work. 


V. V. RESULTS OF TEST RUNS AND 
DISCUSSION 

A complete quantum annealing JSP solver designed to 
solve an instance to optimality using our proposed for¬ 
mulation will require the independent solution of several 
embedded instances {Hp}, each corresponding to a dif¬ 
ferent timespan T. Assuming that the embedding time, 
the machine setup time, and the latency between subse¬ 
quent operations can all be neglected, due to their being 
non-fundamental, the running time T of the approach 
for a specific JSP instance reduces to the expected total 
annealing time necessary to find the optimal solution of 
each embedded instance with a specified minimum target 
probability ~ 1. The probability of ending the anneal¬ 
ing cycle in a desired ground state depends, in an essen- 
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tially unknown way, on the embedded Ising Hamiltonian 
spectrum, the relaxation properties of the environment, 
the effect of noise, and the annealing profile. Under¬ 
standing through an ab initio approach what is the best 
computational strategy appears to be a formidable un¬ 
dertaking that would require theoretical breakthroughs 
in the understanding of open-system quantum anneal¬ 
ing [25] Snjj as well as a tailored algorithmic analysis 
that could take advantage of the problem structure that 
the annealer needs to solve. For the time being, and for 
the purposes of this work, it seems much more practical 
to limit these early investigations to the most relevant 
instances, and to lay out empirical procedures that work 
under some general assumptions. More specifically, we 
focus on solving the CSP version of JSP, not the full 
optimization problem, and we therefore only benchmark 
the Hamiltonians with T = T with the D-Wave machine. 
We note however that a full optimization solver can be 
realized by leveraging data analysis of past results on 
parametrized ensembles and by implementing an adap¬ 
tive binary search. Full details can be found in a longer 
version of this work m- 

On the quantum annealer installed at NASA Ames (it 
has 509 working qubits; details are presented in [32]), 
we run hundreds of instances, sampling the ensembles 
JV = MG {3,4,5,6}, 0 6 {0.5,1}, andP p G {[1,1], [0,2]}. 
For each instance, we report results, such as runtimes, 
at the most optimal Jp among those tested, assuming 
the application of an optimized parameter-setting pro¬ 
cedure along the lines of that described in the previous 
section. Figure [4}a displays the total annealing repeti¬ 
tions required to achieve a 99% probability of success 
on the ground state of Hp, with each repetition lasting 
t-A = 20 ps, as a function of the number of qubits in the 
embedded (and pruned) Hamiltonian. We observe an ex¬ 
ponential increase in complexity with increasing Hamil¬ 
tonian size, for both classes of problems studied. This 
likely means that while the problems tested are small, the 
analog optimization procedure intrinsic to the D-Wave 
device’s operation is already subject to the fundamental 
complexity bottlenecks of the JSP. It is, however, pre¬ 
mature to draw conclusions about performance scaling 
of the technology given the current constraints on cali¬ 
bration procedures, annealing time, etc. Many of these 
problems are expected to be either overcome or nearly so 
with the next generation of D-Wave chip, at which point 
more extensive experimentation will be warranted. 

In Figure [4}b, we compare the performance of the 
D-Wave device to two exhaustive classical algorithms in 
order to gain insight on how current quantum anneal¬ 
ing technology compares with paradigmatic classical op¬ 
timization methods. Leaving the performance of approx¬ 
imate solutions for future work, we chose not to explore 
the plethora of possible heuristic methods as we operate 
the D-Wave machine, seeking the global optimum. 

The first algorithm, B, detailed in [35], exploits the dis¬ 
junctive graph representation and a branch-and-bound 
strategy that very effectively combines a branching 



Number of qubits D-Wave machine annealing time (s) 

FIG. 4: a) Number of repetitions required to solve Hp with 
the D-Wave Two with a 99% probability of success (see the 
appendix). The blue points indicate instances with 9=1 
and yellow points correspond to 9 = 0.5 (they are the same 
instances and runtimes used for Figure [3]c). The number of 
qubits on the x-axis represents the qubits used after embed¬ 
ding. b) Correlation plot between classical solvers and the 
D-Wave optimizer. Gray and violet points represent runtimes 
compared with algorithm B, and cyan and red are compared 
to the MS algorithm, respectively, with 9=1 and 9 = 0.5. 
All results presented correspond to the best out of 5 gauges 
selected randomly for every instance. In case the machine 
returns embedding components whose values are discordant, 
we apply a majority voting rule to recover a solution within 
the logical subspace HD ED HU M HU- We observe a devi¬ 
ation of about an order of magnitude on the annealing time 
if we average over 5 gauges instead of picking the best one, 
indicating that there is considerable room for improvement if 
we were to apply more-advanced calibration techniques |32j . 


scheme based on selecting the direction of a single dis¬ 
junctive edge (according to some single-machine con¬ 
straints), and a technique introduced in [36j for fixing 
additional disjunctions (based on a preemptive relax¬ 
ation). It has publicly available code and is considered 
a well-performing complete solver for the small instances 
currently accessible to us, while remaining competitive 
for larger ones even if other classical approaches become 
more favorable m B has been used in [35] to discuss 
the possibility of a phase transition in the JSP, demon¬ 
strating that the random instances with N = M are par¬ 
ticularly hard families of problems, not unlike what is 
observed for the quantum annealing implementation of 
planning problems based on graph vertex coloring [9]. 

The second algorithm, MS, introduced in [20], proposes 
a time-based branching scheme where a decision is made 
at each node to either schedule or delay one of the avail¬ 
able operations at the current time. The authors then 
rely on a series of shaving procedures such as those pro¬ 
posed by m to determine the new bound and whether 
the choice leads to valid schedules. This algorithm is a 
natural comparison with the present quantum annealing 
approach as it solves the decision version of the JSP in a 
very similar fashion to the time-indexed formulation we 
have implemented on the D-Wave machine, and it makes 





use of the same shaving technique that we adapted as 
a pre-processing step for variable pruning. However, we 
should mention that the variable pruning that we im¬ 
plemented to simplify Ht is employed at each node of 
the classical branch and bound algorithm, so the overall 
computational time of MS is usually much more important 
than our one-pass pre-processing step, and in general its 
runtime does not scale polynomially with the problem 
size. 

What is apparent from the correlation plot in 
Figure [4}b is that the D-Wave machine is easily outper¬ 
formed by a classical algorithm run on a modern single¬ 
core processor, and that the problem sizes tested in this 
study are still too small for the asymptotic behavior 
of the classical algorithms to be clearly demonstrated 
and measured. The comparison between the D-Wave 
machine’s solution time for Hj- and the full optimiza¬ 
tion provided by B is confronting two very different algo¬ 
rithms, and shows that B solves all of the full optimiza¬ 
tion problems that have been tested within milliseconds, 
whereas D-Wave’s machine can sometimes take tenths of 
a second (before applying the multiplier factor ~ 2, due 
to the binary search; see the appendix). When larger 
chips become available, however, it will be interesting to 
compare B to a quantum annealing solver for sizes con¬ 
sidered B-intractable due to increasing memory and time 
requirements. 

The comparison with the MS method has a promising 
signature even now, with roughly half of the instances 
being solved by D-Wave’s hardware faster than the MS al¬ 
gorithm (with the caveat that our straightforward imple¬ 
mentation is not fully optimized). Interestingly, the dif¬ 
ferent parametrized ensembles of problems have distinc¬ 
tively different computational complexity characterized 
by well-recognizable average computational time to solu¬ 
tion for MS (i.e., the points are “stacked around horizon¬ 
tal lines” in Figure j^-b), whereas the D-Wave machine’s 
complexity seems to be sensitive mostly to the total qubit 
count (see Figure |4ja) irrespective of the problem class. 
We emphasize again that conclusions on speedup and 
asymptotic advantage still cannot be confirmed until im¬ 
proved hardware with more qubits and less noise becomes 
available for empirical testing. 


VI. VI. CONCLUSIONS 

Although it is probable that the quantum annealing- 
based JSP solver proposed herein will not prove compet¬ 
itive until the arrival of an annealer a few generations 
away, the implementation of a provably tough applica¬ 
tion from top to bottom was missing in the literature, 
and our work has led to noteworthy outcomes we ex¬ 
pect will pave the way for more advanced applications 
of quantum annealing. Whereas part of the attraction 
of quantum annealing is the possibility of applying the 
method irrespective of the structure of the QUBO prob¬ 
lem, we have shown how to design a quantum annealing 


I. Problem / instance parametrization 

I I 

II. Ensemble pre-characterization III. Choice of mapping 
(software) 

IV. Pre-processing 

I 

V. Ensemble pre-characterization 
(hardware) 

I 

VI. Embedding strategy 

\t \ ' 

VII. Running strategy 

I 

VIII. Decoding and analysis 

FIG. 5: I—II) Appropriate choice of benchmarking and clas¬ 
sical simulations is discussed in Section IV. Ill -IV) Mapping 
to QUBO problems is discussed in Sections II and III. V- 
VI) Pre-characterization for parameter setting is described in 
Section VI. VII) Structured run strategies adapted to spe¬ 
cific problems have not to our knowledge been discussed be¬ 
fore. We discuss a prescription in the appendix. VIII) The 
only decoding required in our work is majority voting within 
embedding components to recover error-free logical solutions. 
The time-indexed formulation then provides QUBO problem 
solutions that can straightforwardly be represented as Gantt 
charts of the schedules. 


solver, mindful of many of the peculiarities of the an¬ 
nealing hardware and the problem at hand, for improved 
performance. Figure [5] shows a schematic view of the 
streamlined solving process describing a full JSP opti¬ 
mization solver. The pictured scheme is not intended to 
be complete, for example, the solving framework can ben¬ 
efit from other concepts such as performance tuning tech¬ 
niques |2S and error-correction repetition lattices |39| . 
The use of the decision version of the problem can be 
combined with a properly designed search strategy (the 
simplest being a binary search) in order to be able to seek 
the minimum value of the common deadline of feasible 
schedules. The proposed timespan discrimination fur¬ 
ther provides an adjustable compromise between the full 
optimization and decision formulations of the problems, 
allowing for instant benefits from future improvements 
in precision without the need for a new formulation or 
additional binary variables to implement the makespan 
minimization as a term in the objective function. As 
will be explored further in future work, we found that 
instance pre-characterization performed to fine tune the 
solver parameters can also be used to improve the search 
strategy, and that it constitutes a tool whose use we ex¬ 
pect to become common practice in problems amenable 
to CSP formulations as the ones proposed for the JSP. 
Additionally, we have shown that there is great potential 




9 


in adapting classical algorithms with favorable polyno¬ 
mial scaling as pre-processing techniques to either prune 
variables or reduce the search space. Hybrid approaches 
and metaheuristics are already fruitful areas of research 
and ones that are likely to see promising developments 
with the advent of these new quantum heuristics algo¬ 
rithms. 
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Appendix A: JSP and QUBO formulation 


and in the penalties version we have 
t' -t < Pi. 

The fact that we are allowing equality in the rewards 
version means that h\ will always have more quadratic 
terms than h\ regardless of variable pruning, leading to 
a more connected QUBO graph and therefore a harder 
problem to embed. 

Another important disadvantage is revealed when 
choosing the coefficients rf, a, and /3 in H' r to guar¬ 
antee that no unfeasible solution has energy less than 
—p'(kN — N). This can happen if the penalty that we 
gain from breaking constraints h 2 or h 3 is less than the 
potential reward we get from h' x . The penalty-based for¬ 
mulation simply requires that 77 , a, and /? be larger than 
0. The following lemma summarizes the equivalent con¬ 
dition for the reward-based case. 

Lemma 1. If / 3 /? 7 / > 3 and a > 0, then 

H' r {x) > ~(k N -N), (A3) 

for all x, with equality if and only if x represents a feasible 
schedule. 


In this appendix we expand on the penalty form used 
for the constraints and the alternative reward-base for¬ 
mulation as well as the timespan discrimination scheme. 


1. Penalties versus rewards formulation 


The encoding of constraints as terms in a QUBO prob¬ 
lem can either reward the respecting of these constraints 
or penalize their violation. Although the distinction may 
at first seem artificial, the actual QUBO problem gener¬ 
ated differs and can lead to different performance on an 
imperfect annealer. We present one such alternative for¬ 
mulation where the precedence constraint ([7]) is instead 
encoded as a reward for correct ordering by replacing 
+rjhi (x) with —rj'h'^x), where 


h'i(x) = Y 


E 




k n -i<i<k n 
\ t+pi<t' 


(Al) 


The new Hamiltonian is 


H' t (x) = —rj'h'^x) + ah 2 (x) + (3h 3 (x). (A2) 


The reward attributed to a solution is equal to 7 / times 
the number of satisfied precedence constraints. A feasi¬ 
ble solution, where all constraints are satisfied, will have 
energy equal to —r]'(knr — N). 

The functions hi and h\ differ only by the range of t 
in the rewards version we have 


t' -t> Pi, 


We also found examples that show that these bounds 
on the coefficients are tight. 

The fact that f3/rj' must be greater than or equal to 3 is 
a clear disadvantage because of the issues with precision 
of the current hardware. In Ht we can set all of the 
penalty coefficients (and hence all non-zero couplers) to 
be equal, which is the best possible case from the point 
of view of precision. 


2. Timespan discrimination 

The timespan discrimination that we propose is a spec¬ 
ification to strike a compromise between the information 
obtained from each solver call, and the required preci¬ 
sion for this information to be accurate and obtained 
efficiently. Specifically, we want this extra information 
to help identify the optimal makespan by looking at the 
energy of the solutions. This means breaking the degen¬ 
eracy of the ground states (i.e., the valid solutions) and 
assigning different energy sectors to different makespans. 
To prevent collisions with invalid solutions, these energy 
sectors have to fit within the logical QUBO problem’s 
gap given by A E = min{? 7 , a, /3}. We note that this will 
affect the actual gap (as seen by the hardware) of the 
embedded Ising model. 

Since the binary variables we have defined in the pro¬ 
posed formulation are not sufficient to write a simple ex¬ 
pression for the makespan of a given solution, additional 
auxiliary variables and associated constraints would need 
to be introduced. Instead, a simple way to implement 
this feature in our QUBO formulation is to add a num¬ 
ber of local fields to the binary variables corresponding 
to the last operation of each job, {Ok 1 , ■■ ■, Ok N }- Since 




10 


the makespan depends on the completion time of the last 
operation, the precedence constraint guaranties that the 
makespan of a valid solution will be equal to the com¬ 
pletion time of one of those operations. We can then 
select the local field appropriately as a function of the 
time index t to penalize a fixed number K of the larger 
makespans ranging from T — K +1 to T. Within a sector 
assigned to the time step T, we need to further divide 
AEj- by the maximum number of operations that can 
complete at T to obtain the largest value we can use as 
the local field hj-, i.e., the number of distinct machines 
used by at least one operation in the set of operations 
{O kl , ■ ■ ■, O kN }, denoted by M fina p If K is larger than 
1, we also need to ensure that contributions from various 
sectors can be differentiated. The objective is to assign a 
distinct T-dependent energy range to all valid schedules 
with makespans within [T—K, T], More precisely, we re¬ 
late the local fields for various sectors with the recursive 
relation 


hr -1 


hr 

-Mfinal 


+ e ) 


(A4) 


where e is the minimum logical energy resolvable by the 
annealer. Considering that e is also the minimum local 
field we can use for hr-K+i and that the maximum total 
penalty we can assign through this time-discrimination 
procedure is A E — e, it is easy to see that the energy 
resolution should scale as A E/(M^ nal ). An example of 
the use of local fields for timespan discrimination is shown 
in Figure 1-d of the main text for the case K = 1. 


Appendix B: Computational strategy 

This appendix elaborates on the compilation methods 
and the quantum annealing implementation of a full op¬ 
timization solver based on the decision version and a bi¬ 
nary search as outlined in the main text. 


1. Compilation 


The process of compiling, or embedding , an instance 
for a specific target architecture is a crucial step given 
the locality of the programmable interactions on current 
quantum annealer architectures. During the graph-minor 
topological embedding, each vertex of the problem graph 
is mapped to a subset of connected vertices, or subgraph, 
of the hardware graph. These assignments must be such 
that the the edges in the problem graph have at least one 
corresponding edge between the associated subgraphs in 
the hardware graph. Formally, the classical Hamiltonian 
Eq. |6| is mapped to a quantum annealing Ising Hamil¬ 
tonian on the hardware graph using the set of equa¬ 
tions that follows. The spin operators scq are defined 
by setting s = 1 and using the Pauli matrices to write 
3i = (erf, of, of). The resulting spin variables erf = ±1, 


our qubits, are easily converted to the usual binary vari¬ 
ables Xi = 0,1 with erf = 2xi — 1. The Ising Hamiltonian 
is given by 


H 

Hq 


H e 


H d 


A{t) [H Q + H E ]+B(t)H D , 

(Bl) 

y! Jij a a, 

ij 

Z I | _Z 

‘ ft \(cn,pj)eE{i,j) k ’ 

,(B2) 

fcev(i) 


- E 

tF -Z_Z 

J i,k,k ,a k a k'i 

(B3) 

(k,k')GE(i,i) 



E << 


(B4) 

kev(i) 




where for each logical variable index i we have a cor¬ 
responding ensemble of qubits given by the set of ver¬ 
tices V(i) in the hardware graph with |H(i)| = Nyyy 
Edges between logical variables are denoted E{i,j) and 
edges within the subgraph of V(i) are denoted E(i,i). 
The couplings J,; 7 and local fields hi represent the logi¬ 
cal terms obtained after applying the linear QUBO-Ising 
transformation to Eq. (JgJ) . Jf k k , are embedding param¬ 
eters for vertex V(i) and (k,k') £ E(i,i) (see discus¬ 
sion below on the ferromagnetic coupling). The equation 
above assumes that a local field hi is distributed uni¬ 
formly between the vertices V ( i ) and the coupling J it j 
is attributed to a single randomly selected edge (ai,fij) 
among the available couplers E(i,j), but other distribu¬ 
tions can be chosen. In the actual hardware implementa¬ 
tion we rescale the Hamiltonian by dividing by Jp, which 
is the value assigned to all Jf kk ,, as explained below. 
This is due to the limited range of Jij and hi allowed by 
the machine m- 

Once a valid embeddin g is chosen, the ferromagnetic 
interactions Jf k k , in Eq. (B31 need to be set appropri¬ 
ately. While the purpose of these couplings is to penalize 
states for which (erf) ^ (erf,) for k,k' £ V(i), setting 
them to a large value negatively affects the performance 
of the annealer due to the finite energy resolution of the 
machine (given that all parameters must later be rescaled 
to the actual limited parameter range of the solver) and 
the slowing down of the dynamics of the quantum system 
associated with the introduction of small energy gaps. 
There is guidance from research in physics [TT; [40j and 
mathematics m on which values could represent the op¬ 
timal Jf k k , settings, but for application problems it is 
customary to employ empirical prescriptions based on 
pre-characterization |6J or estimation techniques of per¬ 
formance [25] , 

Despite embedding being a time-consuming classical 
computational procedure, it is usually not considered 
part of the computation and its runtime is not measured 
in determining algorithmic complexity. This is because 
we can assume that for parametrized families of problems 
one could create and make available a portfolio of embed¬ 
dings that are compatible with all instances belonging to 
a given family. The existence of a such a library would re¬ 
duce the computational cost to a simple query in a lookup 
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table, but this could come at the price of the available 
embedding not being fully optimized for the particular 
problem instance. 


2. Quantum annealing optimization solver 


We now detail our approach to solving individual JSP 
instances. We shall assume the instance at hand can 
be identified as belonging to a pre-characterized family 
of instances for minimal computational cost. This can 
involve identifying N, M, and 0, as well as the approxi¬ 
mate distribution of execution times for the operations. 
The pre-characterization is assumed to include a statisti¬ 
cal distribution of optimal makespans as well as the ap¬ 
propriate solver parameters (Jp, optimal annealing time, 
etc.). Using this information, we need to build an ensem¬ 
ble of queries Q = {q} to be submitted to the D-Wave 
quantum annealer to solve a problem H. Each element of 
Q is a triple (tA, R, T) indicating that the query consid¬ 
ers R identical annealings of the embedded Hamiltonian 
Ht for a single annealing time t A . To determine the el¬ 
ements in Q we first make some assumptions, namely, 
i) sufficient statistics: for each query, R is sufficiently 
large to sample appropriately the ensembles defined in 
Eqs. (B7l (B9); ii) generalized adiabaticity: tA is opti¬ 
mized (over the window of available annealing times) for 
the best annealing performance in finding a ground state 
of Ht (i.e., the annealing procedure is such that the to¬ 
tal expected annealing time tAR* required to evolve to a 
ground state is as small as possible compared to the time 
required to evolve to an excited state, with the same 
probability). Both of these conditions can be achieved in 
principle by measuring the appropriate optimal param¬ 
eters R*( q) and t A (q) through extensive test runs over 
random ensembles of instances. However, we note that 
verifying these assumptions experimentally is currently 
beyond the operational limits of the D-Wave Two device 
since the optimal tA for generalized adiabaticity is ex¬ 
pected to be smaller than the minimum programmable 
value [3]. Furthermore, we deemed the considerable ma¬ 
chine time required for such a large-scale study too oner¬ 
ous in the context of this initial foray. Fortunately, the 
first limitation is expected to be lifted with the next gen¬ 
eration of chip, at which point nothing would prevent 
the proper determination of a family-specific choice of 
R* and t A . Given a specified annealing time, the num¬ 
ber of anneals is determined by specifying ro, the target 
probability of success for queries or confidence level, and 
measuring r q , the rate of occurrence of the ground state 
per repetition for the following query: 


R* = 


l0g[l ~ 7-p] 

log[l - r q ]' 


(B5) 


The rate r q depends on the family, T, and the 
other parameters. The minimum observed during pre¬ 
characterization should be used to guarantee the ground 
state is found with at least the specified ro- Formally, 


the estimated time to solution of a problem is then given 

by 


t=E^ 

q6C 


( fog[l ~ H)] \ 
\l°g[l — r q j / 


(B6) 


The total probability of success of solving the problem 
in time T will consequently be Jl q r q . For the results 
presented in this paper, we used R* = 500 000 and t A = 
min(£ J 4) = 20 fas. 

We can define three different sets of qubit configura¬ 
tions that can be returned when the annealer is queried 
with q. 8 is the set of configurations whose energy is 
larger than A E as defined in Section III of the paper. 
These configurations represent invalid schedules. V is the 
set of solutions with zero energy, i.e., the solutions whose 
makespan T is small enough (T < T — K) that they 
have not been discriminated by the procedure described 
in the subsection on timespan discrimination. Finally, 
S is the set of valid solutions that can be discriminated 
(T £ (T — K, Tj). Depending on what the timespan T of 
the problem Hamiltonian Ht and the optimal makespan 
T are, the quantum annealer samples the following con¬ 
figuration space (reporting R samples per query): 


T <T 
T £ ( T-K,T] 

r <t — k 


V, 5 = 0 -> E 0 > AE, (B7) 

V = 0-^E o £ (0,AE], 

(B8) 

8,V,S E o = 0. (B9) 


Condition (B8) is the desired case where the ground 


state of Ht with energy Eq corresponds to a valid sched¬ 
ule with the optimal makespan we seek. The ground 
states corresponding to conditions (B71 and (B91 are 


instead, respectively, invalid schedules and valid sched¬ 
ules whose makespan could correspond to a global mini¬ 
mum (to be determined by subsequent calls). The above- 
described assumption ii) is essential to justify aborting 
the search when case ( |B8| ) is encountered. If R and t A 
are appropriately chosen, the ground state will be prefer¬ 
entially found instead of all other solutions such that one 
can stop annealing reasonably soon (i.e., after a number 
of reads on the order of R*) in the absence of the appear¬ 
ance of a zero-energy solution. We can then declare this 
minimum-energy configuration, found within (0, AE], to 
be the ground state and the associated makespan and 
schedule to be the optimal solution of the optimization 
problem. On the other hand, we note that if K = 0, 
a minimum of two calls are required to solve the prob¬ 
lem to optimality, one to show that no valid solution is 
found for T = E — 1 and one to show that a zero-energy 
configuration is found for T = T. While for cases (B8) 


and (B91 the appearance of an energy less than or equal 


to AE heuristically determines the trigger that stops the 
annealing of Ht, for case (B7), we need to have a pre¬ 


scription, based on pre-characterization, on how long to 
anneal in order to be confident that T < T . While op¬ 
timizing these times is a research program on its own 
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that requires extensive testing, we expect that the char¬ 
acteristic time for achieving condition (B8) when T = T 
will be on the same order of magnitude for this unknown 
runtime. 


3. Search strategy 


The final important component of the computational 
strategy consists in determining the sequence of times- 
pans of the calls (i.e., the ensemble Q). Here we pro¬ 
pose to select the queries based on an optimized binary 
search that makes informed dichotomic decisions based 
on the pre-characterization of the distribution of optimal 
makespans of the parametrized ensembles as described 
in the previous sections. More specifically, the search is 
designed based on the assumption that the JSP instance 
at hand belongs to a family whose makespan distribution 
has a normal form with average makespan (T) and vari¬ 
ance cr 2 . This fitted distribution is the same V p described 
in Figure 2-a of the main text whose tails have been cut 
off at locations corresponding to an instance-dependent 
upper bound T max and strict lower bound T min (see the 
following section on bounds). 

Once the initial T m [ n and T max are set, the binary 
search proceeds as follows. To ensure a logarithmic scal¬ 
ing for the search, we need to take into account the nor¬ 
mal distribution of makespans by attempting to bisect 
the range T max ] such that the probability of find¬ 

ing the optimal makespan on either side is roughly equal. 
In other words, T should be selected by solving the fol¬ 
lowing equation and rounding to the nearest integer: 


erf 


-fmax + \ ~ (T) 

ay/2 


erf 


+ |-(T) 


aV2 

(BIO) 

T — max(l, K) + \- (T)' 


erf/r+jCm +erf f -- — ' 2 

V aV2 J \ crV2 


where erf(s) is the error function. For our current pur¬ 
pose, an inexpensive approximation of the error function 
is sufficient. In most cases this condition means initial¬ 
izing the search at T = (T). We produce a query qo 
for the annealing of Ht■ If no schedule is found (con¬ 
dition (B7l) we simply let T m i n = T. If condition (B9) 
is verified, on the other hand, we update T max to be the 
makespan T of the valid found solution (which is equal to 
T — max(l, K) + 1 in the worst case) for the determina¬ 
tion of the next query The third condition (B8), only 
reachable if K > 0, indicates both that the search can 
stop and the problem has been solved to optimality. The 
search proceeds in this manner by updating the bounds 
and bisecting the new range at each step and stops either 
with condition (B8I or when T = T max = T m i n +1. Figure 
[6ja shows an example of such a binary search in practice. 
The reason for using this guided search is that the average 
number of calls to find the optimal makespan is dramati¬ 
cally reduced with respect to a linear search on the range 


(T m i n , T max ]. For a small variance this optimized search 
is equivalent to a linear search that starts near T = (T). 
A more spread-out distribution, on the other hand, will 
see a clear advantage due to the logarithmic, instead of 
linear, scaling of the search. In Figure [6}b, we compute 
this average number of calls as a function of N, 9, and K 
for N = M instances generated such that an operation’s 
average execution time also scales with N. This last con¬ 
dition ensures that the variance of the makespan grows 
linearly with N as well, ensuring that the logarithmic 
behavior becomes evident for larger instances. For this 
calculation we use the worst case when updating T max 
due to condition (B9) being met. We find that for the 
experimentally testable instances with the D-Wave Two 
device (see Figure 3-b of the main text), the expected 
number of calls to solve the problem is less than three 
(in the absence of pre-characterization it would be twice 
that), while for larger instances the size of Q scales log¬ 
arithmically, as expected. 



15 20 25 30 35 

Optimal makespan 



FIG. 6: a) View of a guided binary search required to iden¬ 
tify the minimum makespan over a distribution. The fitted 
example distribution corresponds to N = M = 8, fitted to a 
Gaussian distribution as described in the main text. We as¬ 
sume K = 1. The first attempt queries H 26 , th e seco nd H 29, 
and the third H 30 (the final call), following Eq. (BIO I. b) Av¬ 
erage number of calls to the quantum annealer required by the 
binary search assuming Eq. ( |B10[ ) (left panel) or assuming a 
uniform distribution of minimum makespans between trivial 
upper and lower bounds. Thick and dashed lines correspond 
to 6 = 1 and 9 = 0.5, respectively, and the numeric values as¬ 
sociated with each color in the figure correspond to different 
values of K. The operations’ execution times are distributed 
uniformly with P v = {0,..., IV/2}. 




























13 


4. JSP bounds 

The described binary search assumes that a lower 
bound T m ; n and an upper bound T max are readily avail¬ 
able. We cover their calculation for the sake of complete¬ 
ness. The simplest lower bounds are the job bound and 
the machine bound. The job bound is calculated by find¬ 
ing the total execution time of each job and keeping the 
largest one of them, put simply 


benchmark for comparison and performance. Classi¬ 
cal algorithms can sometimes be repurposed as useful 
pre-processing techniques as demonstrated with variable 
pruning. We provide a quick review of the classical meth¬ 
ods we use for this work as well as some details on the 
classical solvers to which we compare. 

1. Variable pruning 


max 

ne{l,...,iV} 


k„ 

X 

i=k n —i 


(BH) 


where we use the lexicographic index i for operations 
and where ko = 1. Similarly, we can define the machine 
bound as 


max 

m£{l, M} 


X ft> 

ielm 


(B12) 


where I m is the set of indices of all operations that need 
to run on machine m m . Since these bounds are inex¬ 
pensive to calculate, we can take the larger of the two. 
An even better lower bound can be obtained using the 
iterated Carlier and Pinson (ICP) procedure described in 
the window shaving subsection of Section III of the main 
text. We mentioned that the shaving procedure can show 
that a timespan does not admit a solution if a window 
closes completely. Using shaving for different timespans 
and performing a binary search, we can obtain the ICP 
lower bound in 0(N 2 \og(N)M 2 T max log 2 (T max UriiT) )) , 
where T m ; n and T max are some trivial lower and upper 
bound, respectively, such as the ones described in this 
section. Given that the search assumes a strict bound, 
we need to decrease whichever bound we chose here by 
one. 

As for the upper bound, an excellent choice is provided 
by another classical algorithm developed by Applegate 
and Cook m for some finite computational effort. The 
straightforward alternative is given by the total work of 
the problem 


X ft- ( B13 ) 

ie{i,fcjv} 

The solver’s limitations can also serve to establish prac¬ 
tical bounds for the search. For a given family of prob¬ 
lems, if instances of a specific size can only be embedded 
with some acceptable probability for timespans smaller 
than T^ 1 ™^ ed , the search can be set with this limit, and 
failure to find a solution will result in at which 

point the solver will need to report a failure or switch to 
another classical approach. 


Appendix C: Classical algorithms 

When designing a quantum annealing solver, a sur¬ 
vey of classical methods provides much more than a 


Eliminating superfluous variables can greatly help mit¬ 
igate the constraints on the number of qubits available. 
Several elimination rules are available and we explain be¬ 
low in more detail the procedure we used for our tests. 

The first step in reducing the processing windows is to 
eliminate unneeded variables by considering the prece¬ 
dence constraints between the operations in a job, some¬ 
thing we refer to as simple variable pruning. We define 
Ti as the sum of the execution times of all operations pre¬ 
ceding operation Oi. Similarly, we define q,; as the sum of 
the execution times of all operations following Oi. The 
numbers r* and qi are referred to as the head and tail 
of operation Oi, respectively. An operation cannot start 
before its head and must leave enough time after finish¬ 
ing to fit its tail, so the window of possible start times, 
the processing window, for operation O, is [Vj,T— Pi — qi\. 

If we consider the one-machine subproblems induced 
on each machine separately, we can update the heads 
and tails of each operation and reduce the processing 
windows further. For example, recalling that Ij is the 
set of indices of operations that have to run on machine 
Mj, we suppose that a,b £ Ij are such that 

r a + Pa + Pb + Qb > T. 

Then O a must be run after Of,. This means that we can 
update r a with 


r a = ma,x{r a ,r b +p b }. 

We can apply similar updates to the tails because of the 
symmetry between heads and tails. These updates are 
known in the literature as immediate selections. 

Better updates can be performed by using ascendant 
sets, introduced by Carlier and Pinson in [36| . A subset 
X C Ij is an ascendant set of c £ Ij if c ^ Ij and 

min r a + 'V' p a + min q a > T. 

aGXU{c} a£X 

a£XU{c} 

If X is an ascendant set of c, then we can update r c with 



• 


j 

r c = max < 

r c , max 
X'cx 

min r a + Pj 

aex' J 

L aex 1 J 

I 


Once again, the symmetry implies that similar updates 
can be applied to the tails. 

Carlier and Pinson in m provide an algorithm to 
perform all of the ascendant-set updates on Mj in 
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0(N\og(N )), where N = \Ij\. After these updates have 
been carried out on a per-machine basis, we propagate 
the new heads and tails using the precedence of the op¬ 
eration by setting 


n+ 1 = max {r i+ i, + p t } , 


(Cl) 


= max { q { , q i+ i + Pi+i} , 


(C2) 


for every pair of operations Oi and O j+ i that belong to 
the same job. 

After propagating the updates, we check again if any 
ascendant-set updates can be made, and repeat the cycle 
until no more updates are found. In our tests, we use an 
implementation similar to the one described in m to do 
the ascendant-set updates. 

Algorithm [T] is pseudocode that describes the shav¬ 
ing procedure. Here, the procedure UpdateMachine(j) 
updates heads and tails for machine i in 0(Nlog(N)), 
as described by Carlier and Pinson in |21| . It returns 
True if any updates were made, and False otherwise. 
Propagate Windows () is a procedure that iterates over the 
tasks and checks that Eqs. (Cl) and (C2) are satisfied, 
in O(NM). 


Algorithm 1 Shaving algorithm 
1: procedure icp_ shave 
2 : updated «— True 

3: while updated do 

4: updated «— False 

5: for i £ machines do 

6: updated «— UpdateMachine(i) V updated 

7: if updated then PropagateWindowsQ 


For each repetition of the outermost loop of Algorithm 
|T] we know that there is an update on the windows, which 
means that we have removed at least one variable. Since 
there are at most NMT variables, the loop will run at 
most this many times. The internal for loop runs ex¬ 
actly M times and does work in 0(N\og(N)). Putting 
all of this together, the final complexity of the shaving 
procedure is 0(N 1 2 3 M 2 T\og(N)). 


2. Classical algorithm implementation 


Brucker et al.’s branch and bound method |.'i5| remains 
widely used due to its state-of-the-art status on smaller 
JSP instances and its competitive performance on larger 
ones [23]. Furthermore, the original code is freely avail¬ 
able through ORSEP [22 ■ No attempt was made at opti¬ 
mizing this code and changes were only made to properly 
interface with our own code and time the results. 

Martin and Shmoys’ time-based approach [20] is less 
clearly defined in the sense that no publicly available 
standard code could be found and because a number 
of variants for both the shaving and the branch and 
bound strategy are described in the paper. As covered in 
the section on shaving, we have chosen the 0(nlog(n)) 
variants of heads and tails adjustments, the most effi¬ 
cient choice available. On the other hand, we have re¬ 
stricted our branch and bound implementation to the 
simplest strategy proposed, where each node branches 
between scheduling the next available operation (an op¬ 
eration that was not yet assigned a starting time) imme¬ 
diately or delaying it. Although technically correct, the 
same schedule can sometimes appear in both branches 
because the search is not restricted to active schedules, 
and unwarranted idle times are sometimes considered. 
According to Martin and Shmoys, the search strategy 
can be modified to prevent such occurrences, but these 
changes are only summarily described and we did not 
attempt to implement them. Other branching schemes 
are also proposed which we did not consider for this 
work. One should be careful when surveying the liter¬ 
ature for runtimes of a full-optimization version based 
on this decision-version solver. What is usually reported 
assumes the use of a good upper bound such as the one 
provided by Applegate and Cook [22]. The runtime to 
obtain such bounds must be taken into account as well. 
It would be interesting to benchmark this decision solver 
in combination with our proposed optimized search, but 
this benchmarking we also leave for future work. 

Benchmarking of classical methods was performed on 
an off-the-shelf Intel Core 17-930 processor clocked at 
2.8 GHz. 
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